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Abstract 

The black hole mass function of supermassive black holes describes 
the evolution of the distribution of black hole mass. It is one of the 
primary empirical tools available for mapping the growth of super- 
massive black holes and for constraining theoretical models of their 
evolution. In this review we discuss methods for estimating the black 
hole mass function, including their advantages and disadvantages. We 
also review the results of using these methods for estimating the mass 
function of both active and inactive black holes. In addition, we review 
current theoretical models for the growth of supermassive black holes 
that predict the black hole mass function. We conclude with a discus- 
sion of directions for future research which will lead to improvement 
in both empirical and theoretical determinations of the mass function 
of supermassive black holes. 



1 INTRODUCTION 

Understanding how and when supermassive black holes (SMBHs) grow is 
currently of central importance in extragalactic astronomy. A significant 
amount of empirical work has established correlations between SMBH mass 
and host galaxy spheroidal properties, such as luminosity [93l 11081 I109j . 
stellar velocity dispersion (the Mbh-c^* relationship, e.g., [HI 11161 1162j ). 
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concentration or Sersic index [501 SS] > bulge mass \101\ 11041 165] , and bind- 
ing energy [H [72] . These scaling relationships imply that the evolution of 
spheroidal galaxies and the growth of SMBHs are intricately tied together. 
The currently favored mechanism for linking the growth of SMBHs and their 
hosts is black hole feedback, whereby black holes grow by accreting gas in 
so-called "active" phases, possibly fueled by a major merger of two gas-rich 
galaxies, until feedback energy from the SMBH expels gas and shuts off the 
accretion process [1451 HOl 1118] . Alternatively, it has been suggested that 
the origin of the scaling relationships does not necessarily require SMBH 
feedback, but emerges from the stochastic nature of the hierarchical assem- 
bly of black hole and stellar mass through galaxy mergers. [1281 178]. 

Feedback-driven 'self-regulated' growth of black holes has been able to 
reproduce the local Mbh-o'* relationship in smoothed particle hydrodynam- 
ics simulations |38 1 11591 180j. Moreover, AGN feedback has also been invoked 
as a means of quenching the growth of the most massive galaxies [201 [33] . 
There have been numerous models linking SMBH growth, the quasar phase, 
and galaxy evolution [Ml [HSl [M [Ml [JTSl [261 EH [JSHl [H]. While feed- 
back is likely important for regulating the growth of SMBHs and galaxies, 
the fueling mechanisms that contribute to growing the SMBH are likely di- 
verse. Major-mergers of gas-rich galaxies may fuel quasars at high redshift, 
and grow the most massive SMBHs. However, major-mergers alone do not 
appear to be sufficient to reproduce the number of X-ray faint AGN [106] . 
and accretion of ambient gas via internal galactic processes [291 EH] ; rnay fuel 
these fainter, lower Mbh AGN at lower z. This is supported by the fact that 
many AGN are observed to live in late-type galaxies out to z ~ 1 [611 146] . 
and the X-ray luminosity function of AGN hosted by late-type galaxies sug- 
gests that fueling by minor interactions or internal instabilities represents a 
non-negligible contribution to the accretion history of the Universe [48j . 

The black hole mass function (BHMF) provides a complete census of 
the mass of SMBHs and their evolution. Because of this, the BHMF is 
one of the primary empirical tools available for investigating the growth 
of SMBHs, and for constraining theoretical models for the growth of the 
SMBH population. Because SMBHs and galaxies are thought to be linked 
in their evolution, the BHMF provides insight into the fueling mechanisms 
that dominate black hole growth, and therefore into the role of feedback 
in the evolution of the host galaxy. The BHMF is also an important tool 
in planning future surveys, as it provides an estimate of the distribution of 
SMBH mass expected for the survey. This in turn is important because mass 
is a fundamental quantity of the black hole, and therefore is an important 
observational quantity for empirical studies of black hole accretion physics 
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[1131 H2| 11071 133 11561 l86j . Of course, further improvement to our under- 
standing of black hole accretion physics will further improve our modeling 
and understanding of black hole accretion and feedback, which in turn will 
improve our understanding of black hole-galaxy coevolution. Therefore, the 
BHMF is an important empirical quantity for SMBH studies. 

In this review we discuss the current status of BHMF estimation and 
theoretical modeling. In § [2] we discuss the non-trivial task of estimating 
the BHMF. In §[3] we discuss current estimates of the local SMBH BHMF. 
In § m we discuss BHMF estimates derived by combining the local BHMF 
with the AGN luminosity function via a continuity equation. In § [5] we dis- 
cuss BHMFs estimated for AGN only. In § [6] we review theoretical models 
for SMBH growth that predict the SMBH BHMF. Finally, in §[7] we discuss 
directions for future improvements to the empirical and theoretical studies 
of the BHMF. We note that unlike, say, the luminosity function, the divi- 
sion between 'observational' and 'theoretical' studies is not as clear for the 
BHMF, as some amount of modeling is necessary in order to estimate the 
BHMF from strictly observational quantities. We have attempted to divide 
the studies according to whether the BHMF is constrained empirically, as 
in, say, a formal statistical fitting procedure, or if it is predicted from a theo- 
retical model for SMBH growth. In reality, the line between theoretical and 
empirical studies is blurry, and some procedures which we have considered 
to be empirical may be thought of as theoretical. 

2 Estimating the Black Hole Mass Function 

The black hole mass function, denoted as (j){MBH , z)dMBH , is the number 
of sources per comoving volume V{z) with black hole masses in the range 
Mbh, Mbh + dMBH- The black hole mass function is related to the joint 
probability distribution of Mbh and z, p{Mbh, z), as 



The normalization of the BHMF is N, the total number of SMBHs in the 
observable universe, and is given by the integral of over Mbh and V{z). 

2.1 Complications with Estimating Black Hole Mass Func- 



Similar to luminosity function estimation, the BHMF may be estimated 
from astronomical surveys. However, while there are many well-established 
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methods for estimating luminosity functions, there are two comphcations 
that make BHMF estimation a more difficult problem [87|. The first is the 
issue of incompleteness. Surveys are typically constructed by finding the set 
of objects of interest containing a SMBH that satisfy a flux criteria, e.g., 
all objects brighter than some flux limit. Surveys are not constructed by 
selecting on mass. Because there is a distribution of luminosities at a given 
SMBH mass, whether it be the luminosity of the host galaxy or of the AGN, 
some SMBHs will scatter above the flux limit and some below. This creates 
a selection function which is less sensitive to Mbh , and it is possible that a 
survey may be incomplete in all mass bins. 

The second complication is the large uncertainty in SMBH mass among 
mass estimators. Currently it is not possible to obtain reliable mass esti- 
mators for large numbers of SMBHs through dynamical and modeling of 
the stellar or gaseous components, and thus scaling relationships are em- 
ployed. Masses may be estimated using scaling relationships between Mbh 
and the properties of the host galaxy bulge or the luminosity and the width 
of the broad emission lines for AGN \178\ 1170] . It has also recently been 
suggested that the X-ray variability properties of AGN may also provide 
another scaling relationship for estimating Mbh |1231 11881 l86j , but further 
work is needed for developing this. While these scaling relationships enable 
one to estimate Mbh for large numbers of SMBHs, they also contain a sig- 
nificant intrinsic statistical scatter. Giiltekin, K., et al. |60) find that for 
early type galaxies there is an intrinsic scatter in Mbh of 0.31 it 0.06 dex 
and 0.38 ± 0.09 dex at fixed host galaxy bulge dispersion and luminosity, 
respectively; the amplitude of the scatter is larger for late type galaxies. 
For AGN with broad emission lines, Vestergaard & Peterson |170j estimate 
the scatter in Mbh at fixed luminosity and line width to be ~ 0.4 dex, 
depending on which emission line is used. 

The statistical uncertainty in the mass estimates can have a significant 
effect on the inferred BHMF. The distribution of the mass estimates is the 
convolution of the intrinsic BHMF with the error distribution in the mass 
estimates. In general, it is typically assumed that the error in the mass 
estimates is independent of the actual value of Mbh- This is not the case 
for Mbh estimated through dynamical modeling, however independence be- 
tween Mbh and its error is likely to be a good approximation for Mbh esti- 
mated using scaling relationships. Because scaling relationships are the only 
feasible manner in which to estimate Mbh for a large sample of SMBHs, 
which is necessary for any estimate of the BHMF, we will assume that Mbh 
and its error are independent. Under the assumption of independence be- 
tween the estimated Mbh and its error, the BHMF that would be inferred 
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Figure 1: Illustration of the bias in the estimated BHMF derived from mass 
estimates. Shown is the true mass function (thick solid black line) for a 
simluated sample, and the mass function derived from the mass estimates 
when the statistical error in the mass estimates is 0.3 dex (red dashed line), 
0.4 dex (green dot-dashed line), and 0.5 dex (solid blue thin line). The mass 
function estimated from the mass estimates is biased, especially at the high 
Mbh end and for large statistical error. 

directly from the distribution of the mass estimates is broader than the in- 
trinsic BHMF, and is thus biased. Figure [1] illustrates this effect, where 
an intrinsic mass function is compared with the distribution of an unbiased 
mass estimator having a statistical uncertainty of 0.3, 0.4, and 0.5 dex, re- 
spectively. As can be seen, the distribution of mass estimates is significantly 
different than the intrinsic mass function. In particular, the distribution of 
the mass estimates falls off more slowly with increasing Mbh, and overpre- 
dicts the number of SMBHs at the high Mbh end of the mass function. 
The bias is worse when the dispersion in the scatter in the mass estimates 
becomes larger. 
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2.2 Methodology for Estimating the Black Hole Mass Func- 
tion 



In order to estimate the SMBH mass function in an unbiased manner, it 
is necessary to match the mass function with the observed distribution of 
the mass estimates and any additional observational quantities that the 
selection functioE0 depends on. The basic idea is to start with an assumed 
mass function. Then, calculate the distribution of mass estimates implied by 
this mass function. In addition, calculate the distribution of observational 
quantities that one's sample is selected on, say, flux, that is implied by the 
assumed mass function. This step allows one to correct for incompleteness, 
but requires an additional assumption about how to relate the mass function 
to the quantity that one's sample is selected on. Finally, impose the selection 
function for the sample, and compare the predicted observed distributions 
of mass estimates and any other observables (e.g., flux) with the actual 
distributions. If they are not consistent, then the data rule out the assumed 
mass function and relationship between Mbh and the observable quantities. 

We can make the above procedure more quantitative by deriving the 
likelihood function for the SMBH mass function. Kelly, Vestergaard, & 
Fan [87] derived the likelihood function for the mass function when using 
masses estimated from AGN broad emission lines. They used this likelihood 
function for developing a Bayesian approach to estimating the SMBH mass 
function. Although their method was limited to broad line mass estimates, 
it is straightforward to generalize their formalism for any generic mass esti- 
mator. Denote the black hole mass estimate as Mbh- In addition, denote 
as X the set of observables that one uses to select one's sample. In the 
majority of cases this will be flux at one or more wavelengths. Then, the 
likelihood function for the BHMF based on a sample of n SMBHs is 



[s{e, V')]'^"" n / p{MbhAMbh,u z,,XMX,\MBH,i, z,,,p)p{MBH,i,Zi\e) dMBm 



where the BHMF is related to N and the probability distribution of Mbh 
and z via Equation ([T]). Here, is the binomial coefficient, denotes the 
parameters for the BHMF, ifj denotes the parameters for the distribution in 
X at fixed Mbh and z, and s{9,%l)) is the probability of including a SMBH 
in one's sample as a function 9 and ip. Here, we have assumed that the dis- 

^The selection function is the probability of including a source in one's sample as a 
function of its measured quantities. 
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tribution in the mass estimates at fixed Mbh, z, and X,p{Mbh\Mbh-, X), 
is known, although one could include additional free parameters for this as 
well. The binomial coefficient arises from the fact that the number of objects 
included in one's survey follows a binomial distributiorH with 'trials' and 
probability of success ■0)- The probability of including a SMBH in one's 
survey as a function of the BHMF, s{9,ip), is calculated from the survey 
selection function z) as 



JQ 



p{X\MBH,Z,^)p{MBH.z\e) dMBH dz 



dX. 



(3) 

It is up to the researcher to choose the particular parameteric form for 
the SMBH mass function, the distribution in the mass estimates at fixed 
Mbh,z, and X, and the distribution of the observable that the sample is 
selected on (e.g., flux) at fixed Mbh and z. Typical choices are log-normal 
distributions, Schechter functions, and mixtures of log-normal distributions. 
Once one has done this, one can use Equation ([2]) to compute a maximum- 
likelihood estimate for the BHMF, or perform Bayesian inference. 

An alternative form of estimating the BHMF can be used when the mass 
estimates are derived from an observational quantity, Y, and the intrinsic 
distribution of Y is known. This is commonly used to estimate the local 
mass function using host-galaxy scaling relationships [105j . In this case, the 
mass function is 



p{Mbh\Y)(I){Y) dY, (4) 

- min 



where (l){Y) is the comoving number density of SMBHs as a function of the 
quantity Y. When both p{A'Ibh\Y) and (j){Y) are known then the BHMF 
follows directly from Equation (jl]). As an example, if the mass function 
is derived from the scaling between Mbh and host galaxy spheroidal lu- 
minosity, Lsph, then Y = Lgph, p{MBH\y) is the Mbh-LspH relationship, 
and (piY) is the luminosity function of stellar bulges hosting SMBHs. As 
with BHMFs determined from a mass estimator, improper treatment of the 
intrinsic scatter in Mbh at fixed Y will lead to a biased estimate of the 
BHMF. However, when calculating the BHMF from Equation (j4]), ignoring 



^Often in the luminosity function literature the likelihood is assumed to be a Poisson 
distribuiton. A poisson distribution is an approximation to the binomial distribution 
when N ^ oo and s{9,'4^) — >■ 0, so Equation converges to the Poisson distribution 
when n <^ N. See [85] for further details. 
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the intrinsic scatter results in an estimated BHMF that is too narrow, un- 
derpredicting the number of SMBHs at the high mass end of 4>{Mbh)- This 
is opposite to the case when one estimates the BHMF directly from mass 
estimates. 

3 Black Hole Mass Functions Derived from Host 
Galaxy Scaling Relationships 

The observed scaling between Mbh and the properties of the SMBH host 
galaxy bulge have motivated several groups to estimate the local BHMF 
p6| fWl El ]M\ ESI Ml Sa [Ml [ESI [na [HZl [ITT], with decreasing 
statistical uncertainties. These estimates of the local BHMF have formed 
the basis for many studies which have attempted to map black hole growth 
by comparing with the AGN luminosity function; this is further discussed 
in § [H Typically, the local BHMF is estimated using the local Mbh-o'* 
relationship or the local MBH~Lsph relationship, combined with the local 
number density of galaxies as a function of stellar velocity dispersion or 
bulge luminosity. 

The scaling relationships between Mbh and host galaxy properties are 
only determined for the local universe, and thus most authors have lim- 
ited their determination of the BHMF based on them to the local BHMF. 
An exception is Tamura, Ohta, & Ueda |161j . who estimated the BHMF 
out to z ~ 1 assuming that evolution in the MBH~Lsph relation is driven 
only by passive evolution in Lsph- Evolution in the scaling relationships is 
currently an area of intense study, with most groups finding evidence that 
the normalization of the scaling relationships increases towards higher z 
[T631 [T29l [TM [THOl [TT5l [9], at least for active SMBHs. However, there are 
still concerns regarding potential biases due to selection effects [97j , but see 
Treu et al. |163] and Bennert et al. [9] for procedures aimed at modeling 
and correcting for selection. There may also be biases due to extrapolating 
the AGN mass estimates derived from the broad emission lines to luminous 
quasars at high z [151]. As such, the uncertainties on the quantitative form 
of the evolution in the scaling relationships and their scatter are currently 
large, limiting their use for determining the BHMF outside of the local uni- 
verse. 

When the Mbh^ct^, relationship is used to estimate the local BHMF, 
it is common to use the velocity dispersion distribution derived from the 
SDSS by Sheth et al. |153j . with an additional component representing the 
brightest cluster galaxies [9S]. Sheth et al. [153] estimate the velocity dis- 
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Figure 2: Local BHMF. The shaded region defines the spread in estimates 
obtained using the Mbh-ct*, MsH^Lsph-, and M^n-Mstar relationships, as 
compiled by Shankar et al. |147j . Based on this estimate the local universe 
is dominated by SMBHs with Mbh < W'^Mq. 

persion distribution for late type galaxies by using the Tully-Fisher relation 
to convert the luminosity function of late type galaxies to a circular veloc- 
ity distribution, and then set a = vd \/2. When the MBH~Lsph relation is 
used it is typical to estimate the distribution of Lgph seperately for early 
and late type galaxies by converting their respective luminosity functions to 
spheroidal luminosity functions using an assumed ratio of bulge luminosity 
to total luminosity. From this it has been inferred that the local BHMF 
is dominated by early type galaxies at Mbh ^ 4 x 10^ |185j . Shankar, 
Weinberg, & Miralda-Escude |147j present a compilation of recently deter- 
mined local BHMFs based on a variety of methods, scaling relations used, 
and data sets used. In Figure[2]we show the range of local BHMFs estimated 
from the Mbr-c*, MBH~Lsph-, and MB_f/-Mstar relationships, as presented 
in Shankar et al. [I37J. In general, estimates of the local mass density of 
SMBHs span the range pbh = (3.2-5.4) x IO^Mq Mpc'^ for h = 0.7 [T47] . 
While the procedure for estimating the local BHMF is, in theory, straight- 
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forward, a number of significant systematics remain. First, there is the 
observational difficulty that most BHMFs derived from the Mbh~<7* re- 
lationship are based on SDSS spectra. Unfortunately, the SDSS velocity 
dispersions are based on a fixed aperture, and thus the size of the aperture 
relative to the bulge varies with the apparent size of the galaxy and its incli- 
nation. In addition, the spectral resolution of SDSS spectra is ~ 100 kms~^, 
making it difficult to reliably measure for SMBHs with Mbh < lO^M©. 
Another concern is that the local BHMF is derived by assuming that the 
Mbh-o'* or MsH-Lsph relations are single power-laws with a constant scat- 
ter in Mbh at fixed fi* or Lgp^. However, recent work has shown these 
assumptions to be incorrect. For one, the Mbh^ct* and MBH~Lsph rela- 
tions diverge at the high Mbh end, which Lauer et al. [96] suggest implies 
that the Mbh-o'* relation is not a single power-law. This divergence creates 
an inconsistency in the BHMFs derived from these two scaling relationships 
|96 1ll65j . Similarly, the a-L relationships for the SDSS and dynamical Mbh 
SMBH samples are inconsistent, suggesting a possible selection bias in the 
estimated BHMFs |187l [T3] . The scatter in the Mbh-c* relation is larger 
for spirals |60j . and appears to increase at low Mbh such that most SMBHs 
lie below the Mbh-o'* relation, [571 e.g.]. Several authors have found dif- 
ferences in the slope and scatter of the scaling relations for pseudobulges 
[771 ESI EH ; however, it is unclear that this result is due to differences in the 
perceived bulge velocity dispersions for bulges as compared to pseudobulges 
or due to different scaling relationships. Kormendy, Bender, & Cornell [92] 
argue that Mbh does not correlate with galaxy disks, and only correlates 
weakly, if at all, with pseudobulges. Although there is still much that we do 
not understand about the Mbh and host galaxy scaling relationship, these 
recent results suggest that the scaling relationships are not a single power- 
law with constant intrinsic dispersion in Mbh, representing a significant 
source of systematic uncertainty in the estimated local BHMF, especially at 
the low-mass end. 

4 Black Hole Mass Functions Derived from the 
Local Mass Function and the AGN Luminosity 
Function 

By employing the argument of Soltan |157| . numerous studies have at- 
tempted to estimate the BHMF at a variety of redshifts by comparing the 
accreted mass distribution implied by the quasar luminosity function with 
the local BHMF jWl [Ml [Ml [IIIl [lEl El EH]. These methods em- 



10 



ploy a continuity equation describing the evolution of the number density of 
SMBHs [271 [T55]: 

d(t>M{MBH,t) d<t)M{MBH,t){M{MBH,t)) 

dt + = r^n.erge{MBH,t). (5) 

Here, {M{Mbh ,t)) is the average growth rate of SMBHs as a function of 
Mbh and cosmic age, t{z), and rimergeiMBHit) is the rate at which the num- 
ber density of SMBHs changes due to mergers of black holes, or ejections 
of black holes from their host galaxies due to gravitational recoil. Techni- 
cally, rimer geiMBH,t) Can also include a contribution from SMBHs which 
are created, but this has not been thought to occur over the redshift range 
in which Equation ([5]) is typically applied, i.e., ^ 5. Because the merger 
rate of black holes is currently unknown, many studies that have employed 

Equation set nmergei^BH it) = 0. 

Under the assumption that SMBHs grow during phases of AGN activity, 
AGN demographics in combination with the local BHMF may be used to 
compute (pMiMBH^t). This is because the AGN luminosity function maps 
the accretion history onto SMBHs, and the local BHMF acts as a boundary 
condition on Equation ([5]); it is also possible in principle to include the 
BHMF for AGN, which provides more information. Studies that have used 
Equation ([5]) to estimate the BHMF generally fall into two categories: those 
that assume an AGN lightcurve, and those that employ the BHMF of AGN. 
We discuss each of these seperately. 

4.1 Methods that Assume an AGN Lightcurve 

Most authors employing Equation ([5]) have assumed a parameteric form for 
{M(MBH,t)). The accretion rate is related to the bolometric luminosity 
output of the accretion flow onto the SMBH as L = erMaccC^, where is 
the radiative efficiency of the accretion flow. Mace is the accretion rate of 
matter onto the SMBH, and c is the speed of light. The growth rate of the 
SMBH is M = (1 — €r)Macc, due to the fact that a fraction €r of accreted 
mass is radiated away as energy. Making this substitution, the continuity 
equation becomes 

d(t)M{MBH,t) 1 - er d^M {Mbh , t) {L{Mbh ,t)) ^ 

dt ^ erC^ dMBH ^ ' 

where we have ignored mergers of SMBHs. Equation ([6]) shows that it is pos- 
sible to calculated the BHMF at an time given the local BHMF, an assumed 
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average accretion flow lightcurve as a function of Mbh-, {L{MBH,t)), and 
an assumed radiative efficiency. Because (pMiMsH, z) and {L{MBH-,t{z))) 
imply a luminosity function, the local BHMF and AGN luminosity function 
can be used to place constraints on and {L{Mbh it)) ■ This means that, in 
practice, one also has to assume a bolometeric correction, which itself likely 
depends on both black hole mass [H] and L/Lem [16611167] . In addition, an 
estimate of {L{Mbh ,t)) also enables one to estimate the lifetime and duty 
cycle of AGN activity, modulo some luminosity-dependent definition of an 
AGN; note that the AGN duty cycle defines the fraction of SMBHs that are 
'active' at a given Mbh and z. 

A variety of lightcurve models have been used when employing Equa- 
tion ([6]) to reconstruct the evolution of the BHMF. The simplest model 
is that where SMBHs spend a fraction of their time radiating at a con- 
stant Eddington ratio, and spend the remainder of their time in quies- 
cence. The free parameters in this model are the Eddington ratio, AGN 
lifetime or duty cycle, and radiative efficiency. This model has been used 
by [IMl [Ml [Ml [HSl fTTTpl to study the build-up of the local black hole 
mass function, although |147j also considered models where the average ac- 
cretion rate relative to Eddington falls off toward lower z and higher Mbh- 
Raimundo & Fabian [133j employed a variation on the constant L/LEdd 
models, assuming three different populations of AGN with their own Ed- 
dington ratio: a population of obscured low L/LEdd AGN, a population of 
obscured AGN with higher L/LEdd-, and a population of unobscured AGN. 
Yu & Lu |185j modeled the quasar lightcurve as radiating at the Eddington 
limit for a period of time, and then transitioning into a power-law decay. 
Cao [23] also modeled the quasar lightcurve as undergoing a power-law decay. 
Lightcurves undergoing a power-law decay arise from self-regulation models, 
and describe the evolution of the lightcurve after black hole feedback un- 
binds the accreting gas, therefore quenching its fuel supply. The power-law 
decay occurs either as a result of evolution of a blast-wave [681 [69j or from 
viscous evolution of the accretion disk [1861 [89] . 

In Figure[3|we compare the BHMF calculated by Shankar et al. [147] with 
that calculated by Cao [23]. For Shankar et al. [147] we show their reference 
model, which assumes a radiative efficiency of = 0.065, an accretion 
rate relative to Eddington of M/MEdd = 0.6, and that half of all SMBHs 
are active at z = 6. We show the model from Cao [23] which assumes a 
radiative efficiency of = 0.11 and a quasar lifetime of 7.5 x 10^ yr, as it 

^Technically, [136] assumed that the Eddington ratio was a weakly increasing function 
of luminosity. 
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Figure 3: Comparison of two recently estimated BHMFs, calculated by 
Shankar et al. [147] (red dashed line) with that calculated by Cao [23j (solid 
black line). Both BHMF were estimated by assuming a quasar lightcurve, 
where Shankar et al. |147j used a step function model while Cao [23] as- 
sumed a power-law decay. Despite the two different models, the BHMFs are 
similar at a variety of redshifts, except at possibly the high mass end. 



better matches the Shankar et al. [147j estimates. The two estimates of the 
BHMF agree fairly well, despite the different quasar lightcurve models. 

In general, most of the studies that have used Equation ([6]) in combina- 
tion with an assumed quasar lightcurve have concluded the following: 

• Most SMBH growth occurs in periods when the quasar is radiating 
near the Eddington limit. 

• Most, if not all, of the local black hole mass function can be explained 
as the relic of previous AGN activity, implying that mergers of SMBHs 
are not important for building up the local mass function. 

• SMBH growth is anti-hierarchical, with the most massive black holes 
growing first. This has also been termed 'down-sizing' of active SMBHs. 
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• The lifetime of AGN activity is ~ a few x 10 yr. 

• Most SMBHs have non-zero spin, as imphed by inferred radiative effi- 
ciencies of er > 0.06. 

However, while Equation ([6]) has proven to be an important tool for study- 
ing SMBH growth, and estimating the black hole mass function, it must be 
kept in mind that use of Equation ^ often entails some strong assumptions. 
These methods rely on the assumed form of the quasar lightcurve, distri- 
bution of radiative efficiencies, and bolometeric corrections, all of which are 
subject to considerable uncertainty. Moreover, in general these methods 
also rely on an estimate of the local black hole mass function, which, as dis- 
cussed in § [3l is itself subject to considerable uncertainty. Indeed, there is 
a strong degeneracy between the estimated radiative efficiency of accretion 
and the normalization of the local BHMF, and therefore the uncertainty in 
€r is linearly proportional to that in the normalization (or integral) of the 
local BHMF. All of these issues have the potential to introduce systematic 
error into methods based on Equation and further work is needed in 
reducing these systematics. 

4.2 Methods that Include the Distribution of Active Super- 
massive Black Holes 

An alternative to the methods described in § 14.11 is to estimate the average 
value of the accretion rate onto SMBHs directly from the observational data. 
This avoids the issue of assuming a form for the quasar lightcurve, as instead 
{L{Mbh it)) is derived directly from the estimated distribution of L/LEdd- 
Techniques based on this approach require a means of linking the mass 
function of active SMBHs to observational quantities, which is done via 
scaling relationships. This was the approach of Merloni [lllj and Merloni 
&: Heinz [112] . who employed the black hole 'fundamental plane' (BHFP) 

[ml ED. 

The BHFP is a scaling relationship between Mbh, radio luminosity, 
and X-ray luminosity, that exists for low-accretion rate black holes (i.e., 
M / MEdd i$ 0.01), extending from galactic black holes to supermassive ones. 
The BHFP likely reflects the connection between Mbh and the conversion 
of the accretion flow into radiative energy and jet power. It, in principle, 
enables one to connect the radio and X-ray luminosity functions to the mass 
function of active SMBHs. Having obtained a distribution of Mbh and X- 
ray luminosity at a given redshift for the active SMBH population, Merloni 
[lllj and Merloni &: Heinz [112j then convert this to a joint distribution 
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Figure 4: Figure 5 from Merloni &: Heinz [112J, showing the redshift evolu- 
tion of their estimated BHMF. The dashed line is the local BHMF and the 
shaded regions reflect the uncertainty in the BHMF that is due to uncertain- 
ties in the AGN luminosity function. The high mass end of their estimated 
BHMF is built up faster than the low mass end, a phenomenon that has 
been called 'downsizing'. 

of Mbh and Mace assuming a conversion from X-ray luminosity to Mace 
which depends on the Eddington ratio. The joint distribution of Mbh and 
Mace at a given redshift for active SMBHs therefore enables calculation 
of the average growth rate {M{MBH,tiz))), which can then be combined 
with the continuity equation to calculate the black hole mass function at 
the next redshift. Their estimated BHMF is shown in Figure [H which is 
a recreation of their Figure 5. Similar to methods based on assuming a 
quasar lightcurve, Merloni Sz Heinz |112| concluded that SMBHs grow anti- 
hierarchically; however, in contrast to the lightcurve methods, Merloni & 
Heinz [112j concluded that most SMBHs have low spin as inferred from their 
derived radiative efficiency. In addition, Merloni & Heinz [112j concluded 
that the distribution of SMBH accretion rates is broad, and that most SMBH 
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growth occurs during a radiatively efficient accretion mode. 

The method of estimating the BHMF from the BHFP developed by Mer- 
loni [111] and Merloni & Heinz |112j has the advantage that it derives the 
distribution of accretion rates empirically. However, there are also disad- 
vantages to this approach. The uncertainties regarding the bolometeric cor- 
rection, estimation of the local BHMF, and radiative efficiency also apply to 
the BHFP method as well. Moreover, as discussed in Merloni &: Heinz | 112j . 
the BHFP is only defined for low-accretion rate objects, i.e., objects with 
L/LEdd < 10"^. Merloni & Heinz [Il2] extrapolate the BHFP to higher 
accretion rates, after rescaling the normalization to ensure that the radio 
luminosity is weak for AGN in the radiatively efficient mode (those objects 
with L/LEdd ^ 10^^ and lacking a jet). Unfortunately, the AGN in the 
radiatively efficient mode make a significant contribution to the X-ray lumi- 
nosity function, from which {M{MBH,t) is derived. Moreover, most stud- 
ies, including those based on the BHFP have concluded that most SMBH 
growth occurs at L/LEdd ^ 10~^, which corresponds to the radiatively effi- 
cient mode. Because the radiatively efficient mode also corresponds to the 
regime of largest systematic uncertainty for the BHFP, there is the poten- 
tial for significant systematic error in estimating the BHMF based on the 
BHFP, as well as in estimating the primary mode of SMBH growth. There 
is thus a need for further improvement to our understanding of the scaling 
relationships involving Mbh and the AGN SED. 

5 Black Hole Mass Functions of AGN 

Thus far, we have focused on methods for estimating the mass function of 
all SMBHs. In this section we will describe methods for estimating the 
BHMF for those SMBHs in AGN, and the results that have come from the 
application of these methods. 

5.1 Methods Based on Scaling Relationships Involving the 
Broad Emission Lines 

The steady improvement in reverberation mapping of AGN [1301 [10] has 
revealed a correlation between the luminosity of AGN and the broad line 
region radius |8H 111) . It is therefore possible, in principle, to obtain an 
estimate of Mbh for broad line AGN (BLAGN) by combining a luminosity- 
based estimate of the broad line region size with an estimate of the velocity 
dispersion of the broad line region gas obtained from the width of the broad 
emission lines jl78j . These virial mass estimates are then calibrated to the 
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estimates of Mbh obtained from reverberation mapping, which themselves 
are caUbrated to be consistent with the local Mbh~<^* relationship [125^ 1181] . 
Currently, calibrations exist for Ha [53], H/3 [ITOl e.g.], Mg II [TTUl [T69l 
1149] . and C IV [170j . The statistical scatter in the virial mass estimates is 
currently estimated to be ~ 0.4 dex |17U] . although there are indications 
that the scatter may be smaller, at least for the most luminous quasars 
[901 [1501 fT60l [881 IIS]- Moreover, it should be noted that the calibration 
for Mg II is obtained by enforcing consistency in the mean values of the 
Mg II mass estimator and the H/3 and C IV ones, and therefore there is 
currently no direct estimate of the statistical scatter in Mg Il-based virial 
mass estimates. In contrast, the amplitudes of the statistical scatter for H/3 
and C IV are estimated by comparing mass estimates derived from these lines 
with the masses derived from reverberation mapping |170 l. Although there 
is currently very little reverberation mapping data for C IV, the estimate of 
the dispersion in the C IV-based mass estimates should not be biased so long 
as the masses based on reverberation mapping are reliable estimates of the 
true Mbh-, regardless of which emission line was used in the reverberation 
mapping campaign. 

Early estimates of the mass function of SMBHs in BLAGN were obtained 
by binning up the virial mass estimates and applying a 'i/Vmax correction 
[1791 [Ml 11681 1169] . a technique borrowed from luminosity function estima- 
tion. Greene & Ho [54] estimated the local BHMF for BLAGN from the 
SDSS DR4, while Vestergaard et al. [M] estimated the BHMF for BLAGN 
over 0.3 < z < 5 using the uniformly-selected quasar sample from the SDSS 
DR3 [132 . Vestergaard & Osmer [^ estimated the BHMF for the bright- 
est BLAGN using objects from a variety of surveys, as their sample was 
designed to complement the uniformly-selected SDSS DR3 sample. Unfor- 
tunately, as discussed in § 12.11 this method of binning up the mass estimates 
suffers from biases due to the large statistical scatter in the virial mass es- 
timates, and due to the inability of a luminosity-based l/V„iax correction 
to correct for incompleteness in Mbh- Subsequent attempts have further 
improved in their methodology, providing more accurate BHMFs. 

Shen et al. [150] employed a forward-modeling approach where the mass 
function and Eddington ratio distribution were estimated by matching the 
observed distribution of mass estimates and luminosity to that implied by 
the model BHMF and Eddington ratio distribution. Their method accounts 
for incompleteness and the statistical scatter in the mass estimates, but 
lacked statistical rigor in that the matching was done visually. Schulze 
&c Wisotzki [139_j employed a maximum-likelihood technique for estimating 
the local BHMF for BLAGN. Their method corrects for incompleteness in 
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Mbh but does not correct the BHMF for the broadening caused by the 
statistical scatter in the virials mass estimates. Kelly et al. [87] developed 
a Bayesian method that corrects for both the statistical scatter in the mass 
estimates and incompleteness, and used their method to estimate the local 
BHMF of BLAGN from the Bright Quasar Survey [135] • Kelly et al. [SS] 
used the method of [87] to estimate the BHMF of BLAGN at 1 < z < 
4.5 from the mass estimates in the SDSS DR3 quasar sample [168] . The 
BLAGN BHMFs from a variety of studies are compiled in Figure [5l showing 
the evolution of the BHMF from the local universe out io z = 4.5. More 
recently, Shen & Kelly |152j extended the Bayesian method of [57] to include 
a possible luminosity-dependent bias in virial mass estimates derived from 
the emission line FWHM, the existence of which was suggested by Shen 
k Kelly [TM]. Shen & Kelly [152] applied their method to the SDSS DR7 
uniformly-selected quasar sample, independently estimating the BHMF and 
Eddington ratio distribution in different redshifts bins. 

Similar to the methods based on the continuity equation, investigations 
of the BHMF for BLAGN have found evidence for the anti-hierarchical 
growth of SMBHs, i.e., cosmic 'down-sizing' of BLAGN activity. The in- 
ferred Eddington ratio distributions are wide, and the density of SMBHs 
continues to increase toward Eddington ratios which are below the sur- 
vey completeness limit. In addition, Kelly et al. [55] used the BLAGN 
BHMF to estimate the lifetime of broad line quasar activity to be t^i^ ~ 150 
Myr among SMBHs with Mbh ~ 10^ Mq, which is similar to quasar life- 
times inferred from the continuity equation. Kelly et al. [88] also used 
their estimated BHMF to estimate the maximum mass of a SMBH to be 
Mbh ~ 3 X IO^^Mq, which is in agreement with theoretical expectations 

[nsKlis]. 

Mass functions estimated from scaling relationships for BLAGN have the 
advantage that they are derived from estimates of Mbh that are obtained 
for individual sources, providing a more 'direct' estimate of the mass func- 
tion than those based on the continuity equation. However, they have the 
disadvantage that they are only available for a subset of the AGN popula- 
tion, which itself is only a subset of the SMBH population. This compli- 
cates comparison with other SMBH mass functions, as the fraction of AGN 
with broad emission lines is poorly constrained, especially as a function of 
mass. This being said, BHMFs of BLAGN represent a subset of SMBHs 
that are actively growing at the time that they are observed, and, as the 
aforementioned studies have demonstrated, their mass function still contains 
important information on SMBH growth. 

As with all methods of BHMF estimation, the virial mass estimates and 
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Figure 5: Broad-line AGN BHMFs at a variety of redshifts. Shown are the 
local BHMF estimated by Greene &: Ho [S^ (dashed magenta line), Kelly 
et al. [87J (dot-dashed green line), and Schulze &: Wisotzki |139j (solid 
blue line). Also shown are the z > 1 BHMFs estimated by Vestergaard 
et al. |168j (dashed red line) and Kelly et al. [88] (shaded region), where 
the shaded region for the Kelly et al. [88] estimate defines an approximate 
95% confidence region. The BHMFs estimated by Greene & Ho [54] and 
Vestergaard et al. |168j are flux- limited BHMFs, as they did not fully correct 
for incompleteness in Mbh- 



the mass functions derived from them still suffer from systematics. First, 
there is the usual problem of calculating a bolometeric correction, although 
this only affects the estimated Eddington ratio distribution, and not the 
BHMF. Second, there are a few concerns with the virial mass estimates 
which could introduce systematic error; some of these have been discussed 
by Greene & Ho [55]. For one, most of the reverberation mapping data is 
only available for the H/3 line. Because of the limited Mg H data, the Mg 
H scaling relationship is in general not calibrated using objects with black 
hole mass estimates from reverberation mapping. There may be systematic 
effects with luminosity or Eddington ratio when using the FW H M-hased 
scaling relationships [32l I152j , possibly due to a dependence of the broad line 
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region structure on these quantities. Systematic effects on broad line region 
geometry, which can effect the inferred velocity dispersion, are a particular 
concern for C IV, which is thought to arise in an accretion disk wind jl35]. 
Along these lines, unaccounted for radiation pressure on broad line clouds 
may also bias the virial masses, especially among those AGN radiating near 
the Eddington limit, |102j : however, its importance is still debated jl2Hll03"t 
1122] . In addition, the reliability of line width measurements can rapidly 
deteriorate for low S/N data [36j- And finally, the BLAGN virial mass 
estimates are calibrated to the reverberation mapping derived masses, which 
themselves are calibrated to lie on the local Mbh~o'* relationship. Most of 
the AGN that are used to calibrate the reverberation mapping masses to 
the Mbh-(^* relationship have lower masses and are hosted by late type 
galaxies, for which there is evidence that the Mbh~o'* relationship begins 
to break down |57) . Greene et al. [57] argue that the normalization of 
the scaling relationships inferred when limiting the calibration to low mass 
SMBHs hosted in late type galaxies may be about a factor of ~ 1.5 lower 
than that used for the current broad line mass estimates |57j . However, 
dynamical mass estimates exist for two reverberation mapped AGN: NGC 
3227 [SIET] and NGC 4151 [IMlllZ]- In both cases the masses derived from 
dynamical modeling and reverberation mapping agree, so it is unclear if a 
smaller scaling factor is needed for late-type galaxies. These issues show that 
there are still many remaining questions regarding virial masses, highlighting 
the need for further study using high-quality reverberation mapping data. 

5.2 Other Methods for Estimating the Black Hole Mass P\inc- 
tion of AGN 

Before broad line mass estimates, there were two earlier attempts at estimat- 
ing the BHMF for AGN, which we briefly mention here. Siemiginowska & 
Elvis |142j and Hatziminaoglou, Siemiginowska, &: Elvis |66j used a model 
for the AGN lightcurve arising due to thermal-viscous accretion disk in- 
stabilities |141j to calculate the expected distribution of luminosity at a 
given black hole mass. Based on this calculated distribution, they used the 
quasar luminosity function to constrain the quasar black hole mass function. 
Siemiginowska & Elvis [142] found evidence for SMBH downsizing in AGN, 
consistent with later work. 

Franceschini et al. [35j found a tight correlation between Mbh and the 
total radio power observed in a sample of local galaxies. They then used 
their empirical relationship to estimate the local BHMF derived from the 
local radio luminosity function of galaxies. While many of the objects in 
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their sample are not considered AGN in the traditional sense, Fraceschini et 
al. [l5] argue that this correlation is a signature of an advection-dominated 
accretion flow, thought to dominate at low accretion rates relative to Ed- 
dington. Therefore, while these SMBH may not be 'active' in the quasar 
sense, the determination of their mass function relies on radio emission from 
the SMBH accretion flow, so this method may still be considered a method 
for estimating the BHMF for active SMBHs. Franceschini et al. [35] com- 
pared their BHMF to models of AGN activity and found that it was incon- 
sistent with AGN activity being continuous and long-lived, but consistent 
with AGN activity being transient and possibly recurrent. 

6 Theoretical Models for Black Hole Mass Func- 
tions Across Cosmic Time 

There have been numerous theoretical models for the formation and growth 
of supermassive black holes, and coevolution with their host galaxies. Un- 
derstanding this formation, growth, and coevolution is one of the current 
most important outstanding issues in extragalactic astrophysics. Because 
the black hole mass function provides a census of the SMBH population 
and its evolution, it is one of the most fundamental observational quanti- 
ties available for constraining models of SMBH formation and growth. As 
such, many theoretical investigations have predicted a BMHF for compari- 
son with the empirical BHMF. In this section we review some of the models 
for SMBH formation and growth. There have been numerous theoretical 
models for SMBH growth and formation, and it is beyond the scope of this 
primarily empirically-focused review to review all of them; instead, we focus 
on those theoretical models that predict a BHMF. 

6.1 Modeling the Coevolution of SMBHs and Galaxies: Pre- 
dicted BHMFs 

Early models for the coevolution of SMHBs and galaxies linked the growth 
of black holes to the properties of host dark matter halos, with periods of 
SMBH growth occuring in quasar phases initiated by mergers. In general, 
early studies that predicted a BHMF used various perscriptions to relate 
Mbh to the mass of the host halo [Ml ESI EH |28l [39] . More recent models 
for the coevolution of SMBHs and galaxies have incorporated AGN feed- 
back from the SMBH. In addition, the availability of empirical BHMFs have 
enabled modelers to compare their more recent models with observational 
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data. In general, the models are qualitatively in agreement with the empir- 
ical results, in that they are able to match the local BHMF fairly well and 
predict downsizing of SMBHs. However, considering the current systematic 
and statistical uncertainties in the empirical results, it is difficult to place 
rigorous empirical constraints on the models such that certain models may 
be ruled out. Because of this, we simply summarize some of the different 
recent models that have been developed which predict the BHMF. 

Granato et al. [52] developed a model incorporating feedback from AGN 
and supernovae, where the feeding of the SMBH is driven by stellar radiation 
drag on gas. Their predicted local BHMF agrees with that estimated by 
Shankar et al. |146j . Cattaneo et al. p6j used halo merger trees constructed 
from A^-body simulations to track the growth of SMBHs. In their model 
the black hole fueling rate was proportional to the star formation rate of 
the host galaxy burst component and the density of the cold gas in the 
starburst component. Their model predicted SMBH downsizing, with the 
most massive part of the BHMF being built up first, in agreement with the 
subsequent empirical studies. 

Hopkins et al. [70j describe a model for the coevolution of SMBHs and 
galaxies whereby all major mergers of gas-rich galaxies trigger a quasar. 
In this model the final black hole mass is assumed to be on average pro- 
portional to the host spheroidal mass, in agreement with the local scaling 
relationships between SMBHs and their host galaxies. Hopkins et al. [70] 
estimated the merger rate of gas-rich galaxies by combining theoretical con- 
straints of the halo and subhalo mass functions with empirical constraints 
on halo occupation models. Their model also predicts SMBH downsizing, 
and their predicted BHMF matches the local BHMF derived by Marconi et 
al. |105) . Similarly, Shen [148] also assumed that quasars are triggered by 
major mergers of gas-rich galaxies, with the SMBHs growing via accretion 
in these quasar phases. Shen |148j used a halo merger rate based on theoret- 
ical expectations from A^-body simulations, and assumed a universal quasar 
lightcurve shape having an exponential increase followed by a power-law 
decay (see also [185] ). The BHMF predicted by Shen [148] broadly agrees 
with the local one estimated by Shankar et al. [147] and predicts that most 
SMBHs with Mbh > 3 x 10^ M© were in place by z = 1, but only 50% of 
them were assembled by z = 2. 

Most recently, Fanidakis et al. |43|.I44| extended the model of [20!, which 
includes AGN feedback, to also follow the spin distribution of SMBHs. In 
their model SMBHs are fueled through accretion of cold gas from mergers, 
disk instabilities, and cooling flows from hot halos. However, the inclusion 
of SMBH spin enabled them to include different radiative efficiencies, which 
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dictates how much accreted material actually grows the black hole, and to 
provide an improved model for the amount of mechanical feedback imparted 
through an AGN jet, both of which depend on the spin of the black hole. 
Their model predicts that the present-day Universe is dominated by SMBHs 
with Mbh ~ 10^-10^ Mq, and that the BHMF at Mbh > 10^ Mq was largely 
built up at z < 2 due to an increase in both lower accretion rate 'radio-mode' 
growth and mergers of SMBHs. 

Almost all models for the cosmological coevolution of SMBHs and galax- 
ies that predict a BHMF have been of an analytical or semi-analytical nature. 
An exception is the study done by Di Matteo et al. [37], who present the 
results from cosmological hydrodynamic simulations of the ACDM model 
that follow the growth of galaxies and SMBHs, including their feedback 
processes, at z > 1. Direct cosmological simulations such as these should, in 
principle, provide the most accurate results as to the predicted BHMF, and 
for identifying the relevant physical processes that are important in shap- 
ing the BHMF. However, current cosmological hydrodynamic simulations 
suffer from the fact that they cannot resolve processes on physical scales 
corresponding to the SMBH accretion flow. In fact, Di Matteo [37] use a 
gravitational softening length of e = 2.73/i~^ kpc. Instead, Di Matteo [37] 
employ a subresolution model where the accretion onto the SMBH is es- 
timated using a Bondi-Hoyle-Lyttleton parameterization [76l [181 IT] with 
a correction factor to account for the fact that the Bondi radius is not re- 
solved. They assume a radiative feedback energy efficiency of 5% [38], which 
is the only free parameter in their model and required in order to match the 
normalization of the observed local Mbh^ct* relationship. Their calculated 
BHMF at z = 1 matches the local BHMF for Mbh >'2.x 10^ Mq. In addi- 
tion, Di Matteo |37) also find downsizing in their model, in agreement with 
observations, with the high mass end of the BHMF being largely in place 
by z ~ 2. 

In Figure [6] we compile predicted BHMFs from several recent models 
for SMBH formation and growth [TQI [Ml [131 [371 [ES]. In general, the 
models tend to agree to within a factor of a few with regards to the BHMF. 
However, they diverge at Mbh ^ 10^ Mq, where their predicted SMBH 
number densities can differ by over an order of magnitude. 

6.2 Modeling the BHMF of SMBH Seeds 

Recent work has made improvement to models for the BHMF by focusing 
on theoretical modeling of the distribution of seed SMBHs. The discovery of 
quasars at z « 6-7 with Mbh ~ lO^M© [79l 1117] places strong constraints 
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Figure 6: Compilation of BHMFs predicted by several recent models for 
SMBH formation and growth. Shown are the BHMFs predicted by Hopkins 
et al. [To] (solid black line), Shen |148) (dotted blue line), Fanidakis et al. 
[33] (dashed red line), Di Matteo et al. [37] (dot-dashed green line), and 
Volonteri & Begelman [I73j (dashed-dotted-dotted-dotted magenta line). 
In general the number densities predicted by the models agree to within a 
factor of a few, although they diverge at Mbh ^ 10^ M©. Some of these 
authors did not report BHMFs at each redshift shown, so we only show 
those available at each redshift. 

on the formation of SMBH seeds due to the very limited amount of time 
available at that redshift to grow SMBHs. Lodato & Natarajan [99] derive 
the BHMF of SMBH seeds at z ~ 15 that are the result of the collapse of 
pregalactic disks which have not yet been enriched by metals I98j. Black 
holes formed through such a mechanism have masses Mbh ^ W^Mq, while 
black holes which are the remnants of Pop HI stars have Mbh ~ IO^Mq. A 
similar seed black hole formation mechanism is through 'quasi-stars' [8l[7l[5], 
which are also able to produce seed black holes with Mbh ~ 10^ Mq. 

Volonteri, Lodato, & Natarajan [176] describe a model for the growth of 
SMBHs seeded according to the direct collapse model of Lodato & Natara- 
jan [98] with varying formation efficiencies. In addition, they also compared 
the results from this model using SMBHs seeded from Pop III remnants. 
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Volonteri et al. |176] grow SMBHs through major mergers, and force the 
black hole mass after the galaxy merger to scale with the circular veloc- 
ity of the host halo; additional growth is also provided through black hole 
mergers. Their merger trees are based on a Monte Carlo algorithm based 
on the extended Press-Schechter formalism. They find that most signifi- 
cant differences in the local BHMF with respect to black hole formation 
efficiency occur at Mbh < IO^Mq, with the number density of SMBHs with 
Mbh < W^Mq increasing with increasing formation efficiency. Volonteri 

6 Begelman [173] performed a similar analysis as that of Volonteri et al. 
|176j but instead used SMBH seeds formed via quasi-stars. The BHMFs 
calculated by Volonteri & Begelman |173j match those of Merloni & Heinz 
[n2j at the high mass end, at least at z > 2. 

Natarajan &: Volonteri |120j used a growth and seeding model which is 
very similar to that employed by Volonteri et al. [176J. However, they also 
predict the BHMF for broad line quasars, assuming that 20% of quasars 
are unobscured. They compare the BHMF derived from their model at 
1 < z < 4.5 to the BHMF for broad line quasars reported by Kelly et al. 
[88] . and to the BHMF for all SMBHs reported by Merloni & Heinz pT2] . 
Natarajan & Volonteri [120] concluded that seeds from Pop HI stars have 
difficulty reproducing the BLAGN BHMF, especially at high redshift, while 
seeds resulting from the direct collapse of pregalactic disks do better at 
fitting the high mass end of the BLAGN BHMF at z > 2. 

7 Directions for Future Work and Improvement 

Before concluding this review, we present a discussion of possible future 
empirical and theoretical work relevant to BHMF studies. These include: 

• Better characterization of the SMBH-host galaxy scaling re- 
lationships. Currently, the local BHMF is estimated from the distri- 
bution of host galaxy properties assuming that Mbh has a constant 
log-normal scatter about a single-power law scaling relationship. As 
discussed in § [3l recent observations have provided reason to doubt 
this assumption, suggesting that the correlations break down at the 
highest and lowest masses. This will create biases in the BHMF deter- 
mined from the scaling relationships, which in turn will also affect the 
BHMF estimated from the continuity equation. Further direct Mbh 
estimates from dynamical and kinematic modeling should be obtained 
for a variety of galaxy types, especially at the high and low Mbh 
end. The next class of 25+ m telescopes should provide a significantly 
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improved picture of the scaling relationships, thus providing us with 
more accurate estimates of the local BHMF. 

• Improvements to techniques based on the continuity equa- 
tion. Most studies that have invoked the continuity equation to link 
the local BHMF to the AGN luminosity function have assumed a single 
radiative efficiency, which is equivalent to assuming a single black hole 
spin, and a universal AGN lightcurve. Neither of these assumptions 
are likely to be true, and improvements to this type of modeling should 
include a distribution of SMBH spin and AGN lightcurves. In addi- 
tion, we need to better characterize the bolometeric corrections, which 
remain a significant source of systematic uncertainty. The continuity 
equation techniques should also be extended to map the evolution of 
the full joint 3-dimensional distribution of black hole mass, accretion 
rate, and spin. While this will not necessarily have a direct effect 
on estimating the BHMF, it will provide insight into the dominant 
accretion modes experienced by active SMBH and into the dominant 
fueling mechanism for AGN activity, as the spin distribution traces 
the SMBH fueling history H^. 

• Better characterization of scaling relationships involving Mbh 
for AGN. The dominant scaling relationship for estimating Mbh in 
AGN involves the broad emission lines. However, as discussed in § 15.11 
a number of systematic uncertainties remain. In order to reduce these 
systematics, we need to better understand the broad line region geom- 
etry for the different emission lines, and how it scales with luminos- 
ity, which will provide us with a more accurate conversion from line 
width to velocity dispersion. Moreover, accurate characterization of 
the broad line region geometry should remove the need to calibrate the 
scaling relationships to the local Mbh~(^* relationship, which has its 
own set of systematics. Improvements to reverberation mapping cam- 
paigns and modeling [TU [121 11271 1189^ [2T] , as well as increasing the 
number of AGN monitored for reverberation mapping, will be needed 
in order to really understand the systematics involved in the broad 
line mass estimates. 

There is also the need to better characterize the black hole fundamental 
plane. Because the BHFP describes how the emission mechanisms 
responsible for the radio and X-ray flux scale with Mbh, the BHFP 
coefficients depend on these emission mechanisms. However, these 
emission mechanisms depend on the geometry of the accretion state 
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and the existence of a jet, which in turn depends on the accretion 
rate |114] . and therefore the BHFP coefficients wih be different for 
different classes of AGN [9ll [5l ESI]. In particular, the BHFP is 
currently poorly constrained for 'soft-state' galactic black holes and 
radio-quiet AGN. Therefore, to reduce the systematics involved with 
the BHFP it will be necessary to characterize the scaling relationships 
and their scatter for radio-quiet objects. A correlation between the 
radio and X-ray luminosity has been observed for radio-quiet objects 
[95j . implying that a BHFP should also exist for these objects. In 
order to better characterize the BHFP for radio-quiet objects it will 
be necessary to obtain radio detections for a well-defined sample of 
radio-quiet AGN with reliable Mbh estimates and X-ray detections. 

Finally, there has recently been the discovery of scaling relationships 
involving Mbh and the optical [Ml HOO] and X-ray [ml [Ml |86] vari- 
ability properties of AGN. Mass estimators based on these scaling re- 
lationships have not been rigorously developed yet, nor have they seen 
widespread use. However, the existence of these scaling relationships 
implies that the variability properties may offer another avenue for es- 
timating Mbh and BHMFs, which may become increasingly valuable 
in the era of current and future large time-domain surveys, such as 
Pan-Starrs and LSST. 

• Understanding the redshift evolution of scaling relations. From 
the theoretical point of view, it is clear that high-redshift scaling re- 
lations (or the lack thereof) between SMBH and their hosts provide 
unique and powerful constraints to models for AGN feeding and feed- 
back, which cannot be otherwise distinguished (see e.g. Merloni et 
al. |115) and references therein). 

In practical terms, a better understanding of the evolution of scaling 
relations may also be very advantageous for BHMF studies. As we dis- 
cussed above, current technique for BH mass estimation at z^O involve 
un-obscured, broad emission line QSOs. One can argue that, as long 
as we are restricted to just this class of QSOs, we will have to make 
critical assumptions about the properties of a significant part of the 
population to draw conclusions about the full BHMF (if we wanted, 
for example, to compare with "continuity-equation-based" methods). 
On the other hand, large multi-wavelength surveys do and will pro- 
vide a wealth of information on the host galaxies of obscured AGN 
at high redshift, that represent the numerically dominant part of the 
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growing black holes population [22l |25l 11831 H]- Therefore, if we had 
an independent way to put constraints on the nature of the BH-host 
relation for these objects, we could explore the uncharted territory of 
BHMF for obscured AGN (and for the entire population). Such an 
independent information could come, for example, from IR studies of 
broad emission lines which could act as probes of the BH potential 
less affected by obscuration. The first exploratory works pursuing this 
line of research have recently been published, e.g., [2l I137j . 

From the technical point of view, a lot of work of course is needed 
to better understand how reliable these estimators are. Another big 
"technical" challenge of all studies of the evolution of scaling relations, 
is the fact that they require a thorough assessment of the many ob- 
servational biases one encounters in studying high redshift AGN and 
their hosts [140j . 

• Accounting for black hole kicks in theoretical models. Most 
theoretical models for the BHMF do not include recoiling effects caused 
by the merger of two black holes. However, recent theoretical work on 
black hole recoils suggests that black holes can spend a significantly 
large enough amount of time offset from the central region of the host 
galaxy to alter their growth, thereby increasing the scatter about the 
scaling relationships and decreasing the final black hole mass [1721 [TBI 
[Tin fT5]- On the other hand, Volonteri, Giiltekin, & Dotti [T74] and 
Volonteri, Natarajan, Giiltekin [177j find that ejected SMBHs are rare 
at z < 5, especially for massive SMBHs, suggesting that accounting 
for ejected black holes will not make a significant difference in the 
BHMF. Further improvement in our understanding of the effects and 
frequency of black hole recoil will ensure the accurate implementation 
of black hole recoil into models for SMBH growth. 

• Improvements to our understanding of AGN feedback. Most 
current theoretical models for SMBH growth involve AGN feedback, 
and assume a single efficiency for coupling feedback energy to the 
gas; this feedback efficiency is usally treated as a free parameter. An 
improved physical understanding of AGN feedback will improve the- 
oretical models for the BHMF, as the feedback efficiency affects the 
dynamics of the SMBH's fuel supply, and therefore the amount that 
the SMBH accretes as a function of redshift. Recent high-resolution 
hydrodynamic simulations in one dimension \30\ 1154^ [3T] and two di- 
mensions [94:\ I126j have concluded that AGN feedback efficiency in- 
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creases with the Eddington ratio, and that the values are below the 
value of ~ 5% assumed in many current theoretical models for SMBH 
growth. Further improvements to simulations developed for studying 
AGN feedback will lead to a better physical understanding of AGN 
feedback, which will improve theoretical models for SMBH growth 
and the BHMF. 

On the observational side, future X-ray observations should provide 
considerable improvement in our understanding of AGN feedback. X- 
ray spectra are needed in order to determine the toal column density of 
the gas, and thus its kinetic energy flux, which can be compared to the 
energy output of the SMBH. Current X-ray observations from Chandra 
have found evidence that AGN feedback exists in the local universe 
|41j . However, X-ray calorimeters on future X-ray satellites will be 
needed for further improvement as they provide the high throughput 
and spectral resolution needed to measure column densities and ve- 
locities of ionized gas, and consequently the kinetic energy flux, in a 
large sample of AGN across a broad redshift range. 

• Improvements in resolution and sub-resolution modeling for 
direct hydrodynamic simulations. Full hydrodynamic cosmologi- 
cal simulations offer the most promising avenue for providing a phys- 
ically motivated BHMF without free parameters, and for unambigu- 
ously identifying the relevant physical processes in building up the 
BHMF. However, they currently cannot resolve scales relevant to the 
accretion flow onto the SMBH. Numerical codes based on Adaptive 
Mesh Refinement techniques will provide improvement in resolution, 
but it will likely be a while before hydrodynamic cosmological simu- 
lations are able to follow SMBH growth in large cosmological volumes 
while simultaneously resolving the scales relevant for individual black 
holes. In the meantime, further improvement can be made to the sub- 
resolution modeling employed by current hydrodynamic simulations. 

One way of improving current sub-resolution models may be to imple- 
ment the results on AGN feedback based on the type of work described 
in the previous bullet-point. Another improvement is in modifying the 
sub-resolution model for the SMBH accretion flow. Current methods 
assume the Bondi rate combined with a correction factor to account 
for the fact that the temperature and density of the gas are not re- 
solved at the Bondi radius. Not surprisingly, the growth of the SMBH 
is sensitive to how this correction factor is modeled [19]. Moreover, 
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Figure 7: Comparison of empirical estimates of the BHMF (grey shaded 
region) with BHMFs predicted by theoretical models for SMBH formation 
and growth (red shaded region), in both the local universe (left) and at 
z = 2 (right). Also shown are the estimated BHMF of broad line AGN 
only (solid green line), from Schulze & Wisotzki [139| (left) and Kelly et 
al. [88] (right). The empirical estimates of the BHMF are those shown in 
Figures [3] and IU while the theoretical estimates are those shown in Figure 
[6l The shaded regions define the spread in the estimates and models, and 
may be considered to be a crude estimate of their uncertainty. In general, 
the theoretical number densities are consistent with the empirical ones to 
within a factor of a few. 



sub-resolution models based on the Bondi rate neglect the angular 
momentum of the gas, and thus the Bondi rate may not be represen- 
tative of the actual accretion rate onto the SMBH. Hopkins & Quataert 
|71j used high-resolution hydrodynamic simulations to conclude that 
the Bondi rate was a poor estimate of the actual accretion rate onto 
the SMBH, and describe a sub-resolution model which accurately es- 
timated the actual accretion rate in their simulations. In addition. 
Power, Nayakshin, & King |132j suggest an alternative sub-resolution 
model based on an 'accretion particle' to provide a more accurate esti- 
mate of the black hole accretion rate. The implementation of improved 
sub-resolution models for accretion rate and feedback into cosmologi- 
cal hydrodynamic simulations, as well as further improvements to the 
sub-resolution models, will result in more accurate predicted BHMFs, 
allowing a more insightful comparison with empirical BHMFs. 

In this paper we have reviewed current estimates of the SMBH mass 
function, as well as theoretical models for the BHMF. As discussed above, 
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each of the methods for estimating the BHMF has their own set of system- 
atics. In Figure [7] we compare the empirical estimates of the local BHMF 
(defined by the shaded region in Figure [2]) with the BHMFs predicted by 
the theoretical models compiled in Figure [6l In addition, in Figure [7] we also 
compare the empirical BHMFs at z = 2^ as estimated using the lightcurve 
method (shown in Figure [3|) and the black hole fundamental plane (shown 
in Figure H]), with BHMFs predicted by the theoretical models. In both 
Figures we also include the BHMFs for broad line AGN as estimated by 
Schulze & Wisotzki |139j and Kelly et al. [88]. In general, the theoretical 
models are consistent to within a factor of a few with the empirical estimates 
of the BHMF, although there is a large spread in the models and empirical 
estimates at z = 2. Moreover, the estimated number densities of broad line 
AGN are significantly lower than those of all SMBHs, suggesting that only 
a small fraction of SMBHs are active across a broad range in Mbh, except 
for possibly SMBHs at z ~ 2 with Mbh > 10*^ Mq. 

Despite the differences in the methods for estimating the BHMF, and 
the theoretical models, they have lead to a number of common conclusions. 
In particular, the empirical results have presented a picture whereby SMBHs 
grow primarily via accretion in active phases (Eddington ratios L^oi / L ^dd > 
0.01), that quasar activity is a relatively short-lived phenomenon relative 
to the lifetime of the SMBH and host galaxy (i.e., small 'duty cycles' for 
AGN activity), and that SMBH growth is anti- hierarchical with the most 
massive end of the BHMF being built up first. These empirical results are 
qualitatively in agreement with the steadily improving theoretical models. 
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